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Abstract. In this paper, we investigate the use of compactly supported divergence- 
free wavelets for the representation of the Navier-Stokes solution. After reminding the 
theoretical construction of divergence-free wavelet vectors, we present in detail the bases 
and corresponding fast algorithms for 2D and 3D incompressible flows. In order to compute 
the nonlinear term, we propose a new method which provides in practice with the Hodge 
decomposition of any flow: this decomposition enables us to separate the incompressible part 
of the flow from its orthogonal complement, which corresponds to the gradient component of 
the flow. Finally we show numerical tests to validate our approach. 
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1. Introduction 

The prediction of fully-developed turbulent flows represents an extremely challenging field 
of research in scientific computing. The Direct Numerical Simulations (DNS) of turbulence 
requires the integration in time of the nonlinear Navier-Stokes equations, which assumes the 
computation of all scales of motion. However, at large Reynolds number, turbulent flows 
generate increasingly small scales: to be realistic, the discretization in space (and correlatively 
in time) ought to handle a huge number of degrees of freedom, that is in 3D out of the reach 
of available computers. 

Many tentatives have been done or are underway to overcome this problem: one 
can cite the Vortex Methods which are able to generate very thin scales, or Large Eddy 
Simulation (LES) and subgrid-scale techniques which separate the flow into large scales, that 
are explicitly computed, from the small scales, that are parametrized or statistically computed. 



In that context, wavelet bases offer an intermediate decomposition to suitably represent 
the intermittent spatial structure of turbulent flows, with only few degrees of freedom: this 
property is mainly due to the good localization, both in physical and frequency domains, 
of the basis functions. The wavelet decomposition was introduced in the beginning of the 
90s for the analysis of turbulent flows [8, 23, 21]. Wavelet based methods for the resolution 
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of the Navier-Stokes equations appear later [2, 11, 9, 16, 13]. They have also been used 
to define LES-type methods such as the CVS method [10]. Most of the works cited below 
use a Galerkin or a Petrov-Galerkin approach for the 2D vorticity formulation with periodic 
boundary conditions. However, if we want to turn to the 3D case, with non periodic boundary 
conditions, these approaches are no more available. 

An alternative was at the same period firstly considered by K. Urban and after 
investigated by several authors: they proposed to use the divergence-free wavelet bases 
originally designed by Lemarie-Rieusset [19]. Divergence-free wavelet vectors have been 
implemented and used to analyze 2D turbulence flows [1, 15, 29], as well as to compute 
the 2D/3D Stokes solution for the driven cavity problem [26, 27]. Seeing that divergence-free 
wavelets are constructed from standard compactly supported biorthogonal wavelet bases, they 
allow to incorporate boundary conditions in their construction [5, 28]. 

This research direction is of great interest, since divergence-free wavelets provide with 
bases suitable to represent the incompressible Navier-Stokes solution, in two and three 
dimension. Our objective is now to investigate their feasibility and amenability for such 
problem. The first point lies in avoiding the pressure by projecting the equations onto the 
space of divergence-free vectors. This (orthogonal) projection is the well-known Leray 
projector, and it can be computed explicitly in Fourier space, for periodic boundary conditions. 
Unfortunately, as already noted by K. Urban [28], if we want to explicit the Leray operator in 
terms of divergence-free wavelets, since they form biorthogonal bases (and not orthogonal), 
they would not give rise, in a simple way, to the orthogonal projection onto the space of 
divergence-free vectors. 

Nevertheless, we propose in the present paper to investigate the use of divergence-free 
wavelets for the simulation of turbulent flows. Firstly, we remind the basic ingredients of 
the theory of compactly supported divergence-free wavelet vectors, developed by Lemarie- 
Rieusset [19]. In section 3, we present in detail the bases we proposed to implement in space 
dimensions 2 and 3. We will see that the choice of the complement wavelet basis is not 
unique from this construction, and it induces the values of divergence-free coefficients, for 
compressible flows. We discuss the algorithmic implementation of divergence-free wavelet 
coefficients in dimensions 2 and 3, leading to fast algorithms (in O(N) operations where N is 
the number of grid points). 

Section 4 is devoted to the Hodge decomposition of a compressible field, in a wavelet 
formulation: the method we present uses both the biorthogonal projectors on divergence- 
free, and on curl-free wavelets: our method is an iterative procedure, and we will 
experimentally prove that it converges. The last section presents numerical tests to validate our 
approach: nonlinear compression of 2D and 3D incompressible turbulent flows, and Hodge 
decomposition of well chosen examples, such as the nonlinear term of the Navier-Stokes 
equations. 
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2. Theory of divergence-free wavelet bases 

In this section, we review briefly the relevant properties of wavelet bases, that will be 
used for the construction of divergence-free wavelets. Compactly supported divergence-free 
vector wavelets were originally designed by Lemarie-Rieusset, in the context of biorthogonal 
Multiresolution Analyses. We illustrate the construction with the explicit example of splines 
of degree 1 and 2. For more details, we refer to [19, 6, 17, 27]. 

2.1. Multiresolution Analyses (MRA) 

Multiresolution Analyses (MRA) are approximation spaces allowing the construction of 
wavelet bases, introduced by S. Mallat [22]. We begin with the one-dimensional case of 
functions defined on the real line. 



Definition (MRA): A Multiresolution Analysis of L 2 (M.) is a sequence of closed subspaces 
(Vj) jeZ verifying: 

(1) Vj, V; C V j+1 , n ?e z Vj = {0}, M - eZ V i is dense in L W 



(3) (Shift-invariance) There exists a function <f) G Vq such that the family {<p(.—k) ; k G Z} 
form a (Riesz) basis ofVo. 

The function <p in (3) is called a scaling function of the MRA. 

Here, j denotes the level of refinement. By virtue of the dilation invariance property (2) 
above, we can deduce that each space Vj is spanned by {4>j : k ; k G Z} where 4>j,k(x) = 
2 j / 2 (f)(2 j x - k). 

Wavelets appear as bases of complementary spaces Wy. 



where the sum is direct, but not necessarily orthogonal. In this context (called the biorthogonal 
case), the choice of spaces Wj is not unique. The problem of constructing the spaces Wj 
means to find a function ip, called wavelet such that the system — k) ; k G Z} spans 
Wq. Repeated decomposition of Vj yields the multiresolution analysis of Vj with the wavelet 
spaces: 



(2) (Dilation invariance) 





V j+1 = Vi © Wj 



(1) 



Vj = V ®Wt 



which leads, when j — > +oo, to the wavelet decomposition of the whole space: 



+oo 

L 2 (R) = V ($W e 



Divergence-free Wavelets for Navier-Stokes 4 

As a result, we can write any function / £ L 2 (R) in the basis ip jjk ; j > 0, k £ Z}, with 

<Pk = 4>- —k) and ^j, fe = 2 J '/ 2 ^(2 j • -fc): 

/ = Cfc ^ fc + S d i' fc ^ 

k&L j>0 fcez 

Dual bases: Let a pair (0, ip) of scaling function and wavelet, arising from a Multiresolution 
Analysis, be given, then we can associate a unique dual pair (0*, ip*), such that the following 
biorthogonality (in space L 2 ) relations are fulfilled: for all k £ Z and j > 0, 

>=5 kfi , <0|^* fc >=O, < ^#* fc >= 6 jlQ 5 k>0 , <^|0*>=O (3) 

Moreover the dual scaling functions (jf k and the dual wavelets ip* k have the same structure as 
above: <j>* k = <j>*(- - k) and = 2^ 2 ip*(2^ ■ —k). 

Scaling equations and filter design: Since the function -^4>(^) lives in V , there exists a 
sequence (h k ) (also called the low pass filter) verifying: 

By applying the Fourier transform||, (4) rewrites: 

0(20 = m (0fe) 
where m (£) = ^ Sfcez h k e~ %k ^ is the transfer function of the filter (/i fc ). 

Again, because of W-i C Vo, the wavelet satisfies a two-scale equation: 

1 x 
Sf ( 2 



1 V>(f) = X>«(x-fc) (5) 

fcez 



where the coefficients (g^) are called the /n'g/z pass ///ter. Again the Fourier transform of ^ 
expresses with the transfer function n of filter g k as: 

^(2£) = n o (O0(O 
In the same way, the dual functions satisfy scaling equations: 



7# (I) = E fce z K 0*(x - fc), 0*(2£) = m*(£)0*(O 
Following [17, 6], the biorthogonality conditions (3) for the scaling functions imply: 



(6) 



™o(0™o(0 + ™o(f + 7r)^8(C + tt) = 1 
while one can choose, for example, as transfer functions for the associated wavelets: 



no(0 = e" 4 « m£(£ + vr), ng(0 = e" l « m (£ + tt) 
The Fourier transform of a function/ is defined by /(£) = J_°° f(x) e~ lx ^dx 
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which corresponds to: 

9k = hl_ k , g* k = (-I) 1 "* ^ Vfc 

In practice, the filter coefficients hk and ^ are all what is needed to compute the wavelet 
decomposition (2) of a given function. Notice that these filters are finite if and only if the 
functions ij; and are compactly supported. 

Example: symmetric biorthogonal splines of degree 1 

A simple example for spaces Vj are the spaces of continuous functions, which are piecewise 
linear on the intervals [k2~\ (k + 1)2 _J '], for k E Z. In this case we can choose as scaling 
function the hat function <f)(x) = max(0, 1 — \x\). Its transfer function is given by 

M0 = e*(i^) a (7) 
The shortest even dual scaling function associated with is associated with the filter: 

mS(0 = e* (^y^) (2 -cosC) (8) 

The corresponding values of filters (h^) and (hi) are given in table 1. Figure 1 displays the 
scaling functions and their associated wavelets in this case. 




Figure 1. From left to right: the scaling function with its associated symmetric wavelet with 
shortest support, and their duals: the dual scaling function <fi* and the dual wavelet ip*. 



2.2. Decomposition-recomposition algorithm and useful example 

In the context of a biorthogonal Multiresolution Analysis, the wavelet decomposition of a 
given function / G L 2 (M.) (equation (2)), is obtained through the now well known Fast Wavelet 
Transform [22]. We briefly review here the formula that will be useful for the following. 

In practice we begin with an approximation fj of / in some space Vj of the MRA. 
This approximation may be the oblique projection of / onto Vj following the direction 
perpendicular to VJ (also called the biorthogonal projection on Vj); but usually fj means 
an interpolating function of /, associated to nodes {k2~ J ; k £ Z} (see in last section some 
examples of procedure). 
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This approximation fj of function / may be expanded in terms of the scaling basis 

4>j,k = 2 J/2 <f)(2 J x - k) ofVy. 

+00 

fj(x) = 2 3 ' 2 c ^ <K2 J * - k) 

k=— oo 

The wavelet decomposition of fj corresponds to a truncated sum in equation (2) up to the 
level J — 1, meaning that 2~ 3 is the finest scale of approximation: 

j-i 

fj = ^2°k <l>k + ^2Yl d i' k ^i> k ^ 

k&L j=0 kel, 

The wavelet coefficients dj^ are then computed recursively on the level j, from j — J — 1 to 
j — 1 using the decomposition of spaces (1): 

+00 +00 +00 

fj+l = C j+l,k <f>j+l,k — c j,k 4>j,k + dj t k 1pj,k 

k=—oo k=— 00 k=— 00 

(Here denotes the biorthogonal projection of fj onto Vj+i). By biorthogonality (3) one 
has: 

C i,k =< fj+l\<t>j,h > ' d hk =< fj+M],k > 1 C j+l,k =< fj+lW j+ l,k > 

which yield the decomposition formula: for all j = 0, ... J — 1, 

C j,k = Y2i c j+l,£+2k 

dj,k = l2i 9e c j+i,t+2k 

where the filters h\ and g\ arise from the scaling equations (6) of the dual basis functions. 
In the same way, we obtain the reconstruction formula: 



c j+l,k — (hk-2l Cj,i + 9k-u d 



where and are the filters provided by the scaling equations (4, 5) of the primal basis 
functions. 



The computing cost for the whole wavelet decomposition (9) (as well as for the 
recomposition) is about C2 3 operations, where 2 3 is the number of point values f(k2~ 3 ) 
we start with, and C means the length of the filters (h* k for the decomposition, h k for the 
synthesis). 

Example: spline wavelets of degree 1 and 2: Biorthogonal splines provide with wavelet 
bases which are regular, compactly- supported and easy to implement. The scaling functions 
of the associated MRA are standard B-spline bases, and the wavelets are constructed easily, by 
linear combinations of translated B-splines. We focus here on two examples of wavelet bases, 
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00 0i ^1 



Figure 2. Scaling functions and associated even and odd wavelets with shortest support, for 
splines of degree 1 (left) and 2 (right). 

which will be useful for the construction of divergence-free wavelets: splines of degree 1 (V^° 
MRA spaces) and splines of degree 2 (V^ 1 MRA spaces). In both cases we draw the scaling 
functions and the associated wavelets with shortest support (figure 2). In both cases the filters 
are easy to compute (see [17, 6]). Since the support of basis function is very short, the filters 
have a few non zero coefficients. The values of the decomposition and reconstruction filters 
are given on Table 1. 
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Table 1. Decomposition filter {h* k) <jjjp and reconstruction filter (hk, gk) coefficients, 
associated to piecewise linear splines (left) and piecewise quadratic splines (right), verifying 
(1 1) (see hereafter) with shortest supports. 



2.3. Multivariate wavelets 

The above considerations can be extended to multi-D. The simplest way to obtain multivariate 
wavelets is to employ anisotropic or isotropic tensor products of one-dimensional functions. 
To be more precise, we focus on the two dimensional case: let (V^ ) and (Vj) two 
multiresolution analyses of L 2 (R) be given, associated with scaling functions and wavelets 
(0o, ipo) and (0i, ^i); the two dimensional tensor product space Vj ® V) is generated by the 
scaling basis {^o > j,fci(z)0v > k 3 (y); (^1^2) e Z 2 }, where (f) ,j, kl (x) = 2 J / 2 o (2 J x - fc x ) and 
similarly for (f>i,j,k 2 - Then each function fj in Vj ® Vj can be written: 

00 00 

fj(x,y)= E cj MM 2 J <j> Q {2 J x-k l )<j> l {2 J y-k 2 ) (10) 

fel=— 00 k2=— 00 
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The anisotropic 2D wavelets are constructed with tensor products of wavelets at different 
scales { 1 0oj 1 ,fci(a ; )'0i ) j2,A;2 (y)}- For certain choices of the support of the functions may 
be very lengthened. In this case the wavelet decomposition of fj writes: 

(fei,fc 2 )eZ 2 
J-i J-i 

+ E E 2 ° 1+J2)/2 E ^o(2^x - fci) ^(2*y - fc 2 ) 

J1=0j2=0 (fc l! fe 2 )eZ 2 
The anisotropic decomposition is the most easy way to compute a multi-dimensional wavelet 
transform, as it corresponds to apply one-dimensional wavelet decompositions in each 
direction. In the 2D case, this is schematized in figure 3. 



J.k 



im( x ) x <kjM^y) 









































,(*iifca) 



Figure 3. Anisotropic 2D wavelet transform. 

In the isotropic case, the 2D wavelets are obtained through tensor products of 
wavelets and scaling functions or wavelets at the same scale. This produces the following 
decomposition for fj\ 

fj(%,y)= E c ki,k 2 M x - h) My - k 2) 

(fei,fca)eZ a 

+ E ( E d f,kiM ^0J,hi( x ) <Pl,j,k 2 (y) + E 45£*a ^O.i- fc i( X ) V'lj.fca 

J=0 \fci,fe2 fcl,fc2 

+ E d f,kiM ^ojm( x ) i>ijte(y) ) 

k 1 ,k 2 / 
As one can see, this decomposition involves three kinds of wavelets, one following the 
direction x: ^^'°\x,y) = ipo(x) <fii(y), one following the direction y: ty( 0,1 >(x,y) = 
4>q{x) ipi(y) and one in both directions: ^^'^(x, y) = ^o(x) ipi(y). The interest of this 
basis remains in the fact that the size of their support is proportional to 2~ J in each direction, 
i.e. the basis functions are rather isotropic. The principle of the associated decomposition 
algorithm is illustrated by figure 4. 
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Figure 4. Isotropic 2D wavelet transform. 



2.4. Theoretical ground of the divergence-free wavelet vectors 
Let introduce 

Hdiv( R ") = {/ e (^ 2 (K"))7div/ e L 2 (M n ), div/ = 0} 
the space of divergence-free vector functions in K" . 

The construction of compactly divergence-free wavelets in (L 2 (M n )) n , which will correspond 
to Riesz bases of H^ v (E n ), was originally derived by Lemarie-Rieusset [19]: it is based on 
the following proposition, which relates two different multiresolution analyses of L 2 (R) by 
differentiation and integration: 

Proposition: Let (Vj*) a one -dimensional MRA with a derivable scaling function (f> 1 and a 
wavelet ipi, be given. Then, we can build a MRA {Vj) with a scaling function <p and a 
wavelet ipQ verifying: 

Vq = span{(fi (x — k),k E Z} V — span{<pi(x — k), k G Z} 

and 

(j)[(x) = o (x) - M% - 1) = 4 Ms) (H) 

For the refinement polynomials it can be traduced by: 

2 

1 + e ^ 

Equation (11) rewrites for the dual functions 0q, <f>\, and ip\: 

<(x) = cj>l{x + 1) - = -4 ^(z) (12) 

which induces: 
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Example: As an example of functions fulfilling the above proposition, one shall cite the 
piecewise linear spline functions (f> , associated with the piecewise quadratic spline 
functions 0i, tpi introduced in section 2.2, and plotted on figure 2. 

Then, the divergence-free wavelets are explicitly constructed by combining suitable 
tensor products of these functions. For instance, in the 2D case we may have the following 
basis [19]: 



Example: The 2D divergence-free vector scaling function takes the form: 

$ div (xi,x 2 ) 



-0l'(Xi)0i(x 2 ) 



l(Xi) [<p {x 2 ) - <f) { x 2 - 1)] 
■[(/> (Xl) - 0oOl - 1)] 01 O2) 



and the corresponding isotropic vector wavelets are given by the system: 



- 1^1(^1) [0o O2) - 0o(>2 - 1)] 

V'o(£l)0l(>2) 



(0,1) 

div 



{x\,x 2) 



')o(xi) - <j>o{xi - l)]^i(x 2 ) 



^i(xi)ipo(x 2 ) 

-^0(^1)^1(^2) 



It can be easily seen that the dilated and translated functions ^^j v . ^ = 2 J '\P^ v (2 3 'a;i — 

h±, 2 j x 2 —k 2 ) with j, k\,k 2 G Z and 5 G {0, 1} 2 \(0, 0) span the space H(jj v (R 2 ) of divergence- 
free vector functions in R 2 . We represent in figure 5 the three generating functions in the case 
of spline generators of degree 1 and 2 of figure 2. 

More generally, the construction of divergence-free wavelets in R n is carried out by suit- 
able combinations of tensor products of functions O , tp Q , and 0i, ipi, satisfying the above 
proposition (see [19, 27, 28]). These allow to state the following theorem of existence of 
isotropic divergence-free wavelet bases in the general case [19]: 

Theorem: There exist (n — 1) (2 n — 1) vector functions \P^- . G H^ v (M n ) (e G VC n of cardinal 

(2 n — 1), 1 < i < n — 1) compactly supported, such that every vector function u G H^i v (R n ) 
can be expanded in a unique way: 

E d l 

jeZ^n keZ" 

and one has, for a constant C > independent from u = (ui, u 2 , . . . , u n ): 



div,i,i,k div,i,j,k 



EE E 1 div.i.^k 1 

,- 6 Z £ en„k 6 Z" 



1/2 



<C||u| 



L 2 
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Figure 5. Isotropic 2D generating divergence-free wavelets ^^y (left), ^njy (center) and 
^ (right). 



where ||u|| 2 = \\ui\\ 2 L2 . 

These wavelets have already been studied by several authors for 2D analyses of turbulent 
flows [1, 15], and also to solve the Stokes problem in two and three dimensions [26, 27, 28]. 
From now on, we will focus on the 2D and 3D case, and we will see that the expansion of 
compressible flow in terms of divergence-free wavelet bases is not uniquely given. We will 
also present in a practical way, the associated fast algorithms. 

3. Practical implementation of divergence-free wavelets 

In this section, we present in detail the two and three-dimensional divergence-free wavelet 
bases, and we study how to compute the associated fast wavelet transforms. We present 
the constructions in the isotropic multivariate wavelet case (see section 2.3), and also in the 
anisotropic one, which differs somewhat from previous studies. 

In the following, we are given two ID multiresolution analyses (V®) and (V^ 1 ) satisfying 
the proposition of section 2.4. We note (f> , ^ an d (f>i, 4>i their associated (one-dimensional) 
scaling functions and wavelets. 

3.1. Isotropic divergence-free wavelet transforms 

3.1.1. The 2D case The starting point of the construction lies in considering as 2D 
multiresolution analysis of L 2 (R 2 ) 2 the vector space of tensor-products (V^V®) x (V^V^). 
In the isotropic case, the 2D scaling functions $i, $2 and the 2D wavelets \&f , of this MRA 
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$2(^1, x 2 ) 





(t>o(xi)<f>i(x 2 ) 



12 



* 1 1,0) (x 1 ,x 2 ^ 



M^M^) ^(1,0), , 





i>o{xi)(f)i(x 2 ) 



4>i(x 1 )ip {x 2 ) 




(0,1), 



Xi,X 2 ) 





^0(^1)^1(^2) 



tpi(xi)ipo(x 2 ) 




(1,1), 



xi,x 2 ) 







As explained in section 2.3, the functions 

{M' ) ( k (./ ;. /'•-) = 2AI/ ; (2',- ; - fc 1; 2^ 2 - k 2 )} 

with j G Z,k= (ifeijJfea) eZ 2 ,e e {(0, 1), (1, 0), (1, 1)}, z = 1, 2, form a basis of (L 2 ( 
Then, a velocity field u in (L 2 (R 2 )) 2 has the following wavelet decomposition: 

£ E 



u 



d d,0) ^(1,0) 

i,i,k 1 ,/.k 



2,j,k 2j,k i,j,k i,,.k 



I 2 )) 2 . 
(13) 



2,j,k 2,j,k i,j,k i,j,k 2j,k 2,j,k 
By a simple linear change of basis we are able to find out a divergence-free wavelet basis, 
and its complement: 



(1,0) 



(1,0) 



^(1,0) 

div 

^(1,0) 



¥ 2 ^-\[¥^-¥^\x u x 2 -i)} 



(1,0) 



(0,1) 



(0,1) 



^(0,1) 

div 

^(0,1) 



*1 



(0,1) 



(0,1) 



(0,1) 



(zi - l,x 2 )} 



(1,1) 



(1,1) 



div 

n 



(1,1) 



(1,1) 



(1,1) 

2 

(1,1) 



The first functions yield a divergence-free basis (already shown in example of 
section 2.4), and the second ones ty e n are the complement functions corresponding to non 
divergence-free part of the data. Remark that the functions and are not orthogonal. 
Moreover, the choice of the functions ^> £ n is not unique, and this choice has an influence on 
the values of all the coefficients, when applying this transform to a compressible flow. The 
choice of ^ £ we made in this work, leads to very simple formula to obtain the divergence-free 
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Now, the expansion (13) of a given vector function u can be rewritten: 

divj,k div,j,k ^ divj.k div,j,k div,j,k divj.kj 

ieZkeZ 2 

+ YY (d m y + d {0 '\ + d {1 '\ *M) (14) 

Z-^t \ n,j,K n,j,K n,j,K n.j,k n.j,k n,j.k) v 

ieZkeZ 2 

where the new coefficients are directly expressed from the original ones by: 



d (1 ' 0) 
divj,k 



(ddiv) < 



d 



(o,D 
div,j,k 



d (1 'l 
2,j,k 



ij,k 



(d n ) < 



divj,k 



id [1A) - ±d (1>1) 

2 ij,k 2 2,j,k 



n,j.k 



n,j,k 



d 



(1,1) 
n,.7,k 



d (i,o) i d (i,o) 
ij,k 4 2 j,k 



2,j.k 4 |.^k 



2 ij,k 2 2j,k 



4"2,i,fei,fe 2 -l 

±d (0 ' 1} n<n 

4 a lJ,fci-lWPJ 



As one can see, the computation of divergence-free wavelet coefficients d^ . ^ is reduced to 
a very simple linear combination of the standard wavelet coefficients d £ . ^ provided by the 
biorthogonal fast wavelet transform, which makes programming easy. 



Remark: div(\l/l 1 ' ' ) ), drv(^/n' l> ) and div(\E r n ,i ' 1 ) are generating functions of the scalar space 



,(o,ih 



r (i,i) 



V° ® V°. Moreover, for all u, we have 
div u = 

jeZkeZ 2 



Then, the incompressibility condition div u = is equivalent to d e . ^ = 0, for all j, k, e. 



,(o,i) 



d (1 'l div^ 1 ' 1 ) 



For incompressible flows, since the biorthogonal projectors onto the spaces ( Vj ® V®) x 
(V^° <8> V^ 1 ) commute with partial derivatives [19], the coefficients d^. y ^ are uniquely 
determined, by the formula (d^) in equation (18). We present in section 5.2, numerical 
experiments on 2D incompressible turbulent flows. 

Difficulties arise when we want to compute the divergence-free part of a compressible 
flow. Because of the non-orthogonality between the divergence-free basis ) and its 

complement (^), the values of the divergence-free wavelet coefficients depend on the choice 
of the complement basis. We address this problem in the specific section 4, on the Hodge 
decomposition. 



3.1.2. The isotropic 3D case The construction of the 3D divergence-free wavelet bases may 
be obtained in a similar fashion as for the 2D bases. Again, the key -point is to start with a 
vector multiresolution analysis of (L 2 (R 3 )) 3 , of the type 

{V} ® Vf ® Vf) x {Vf ® V} (g> Vf) x (V, ® Vf ® Vj) 
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$ 1 (x 1 ,x 2 ,x 3 ) 



01 (^l) 00 (^2) 00 (£3) 






$5 










000001 



and 21 generating 3D-vector wavelets: 

{tf^ I i = l,2,3, e = (ex, e 2 , e 3 ) with e* = 0, 1 and e ^ (0, 0, 0)} 
For example, we give below the expressions of wavelets q?( 1,0 '°\ ^r( 1 '°'°) j an( j ^r^AO). 



^'°'°\x 1 ,x 2 ,x 3 ) 



^l(Xi)0o(x 2 )0o(x 3 ) 






(1,0,0) 





^00100 





(1,0,0) 






^00001 



»f W) (xi I x J) is) 



^1(^1)^0 (£2)00(^3) 







(1,1,0) 





^0^100 





(1,1,0) 






■00^001 



® ( i' 1 ' l \x 1 ,X2,X 3 ) 



^1(^1)^0(^2)^0(^3) 






(1,1,1) 





^oV^o 




(1,1,1) 








and it goes similarly for (0, 1, 0), (0, 0, 1), (1, 0, 1) and (0, 1, 1). 

Let introduce = {e G {0, l} 3 \ (0,0,0)}. The isotropic wavelet expansion of a given 
function u writes: 

u = E E E K,k *;„.k + <.,k *2., k + ^,,.k *3., k ) as) 

jGZk e Z 3£6 ^ 

Following theorem of section 2.4, there exist 14 kinds of isotropic divergence-free wavelets, 
with arbitrary possible choices concerning the privileged direction of each basis function. In 
the following we do not detail all the expressions we choose, but only some typical ones: 
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-i^lOl)[0o(>2) - 0o(^2 ~ 1)]0 O (^3) 





-j1pl(xi)(f)o(x2)[M x 3) ~ M x 3 ~ 1)] 



i>o(Xl)M x 2)M x 3) 



*<j£f(si,S2,Ss) 



^l(^l)^o(^2)0o(^3) 
-^0(^1)^1(^2)00(^3) 





^div 2 ( X 1, X 2, X 3) 



-§^i(xi)^o(x 2 )[0o(x 3 ) - 0o(^3 - 1)] 

-^o( x l)lpl( x 2)[(j>o( X 3) ~ M x 3 ~ !)] 
^0(^1)^0(^2)01(^3) 



^i^( x U x 2,X 3 ) 



-ipi(x 1 )tp (x 2 )i'o(x 3 ) 







^0(^1)^0(^2)^1(^3) 



^0(^1)^1(^2)^0(^3) 
--00(^1)^0(^2)^1(^3) 

It goes similarly for all basis function: for each e <E il* 3 given, two divergence-free wavelets 
., i = 1, 2 are carried out by linear combination of \&f, m order to satisfy the 

divergence-free condition. The complement wavelet ^ is constructed in order to take care 
of the symmetry. For example, we consider: 



^(1,0,0) 

div,i 



(1,0,0) w lTf (l,0,0) 



^(iA0) () 



< ^ '°) = * 

dlV,2 



(1,0,0) 1 / lTf (l,0,0) 



- ^ 



(1,0,0), 



i,0) 



1)) 



(1,0,0) 



(1,0,0) 
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and similarly for ^^y and , i 



^(1,1,0) 

div,i 



^(1,1,0) 

div,2 



(1,1,0) 



(1,1,0) 



(1,1,0) l/ lTf (1,1,0) 



*1 



(1,1,0), 



1)) 



.1(^(1,1,0), 



(1,1,0), 



1)) 



(1,1,0) 



(1,1,0) 



(1,1,0) 



and similarly for ^^y. and ^[jjy > i 



< ^(i,i,i) 
div,i 



1,2. 



(i,i,i) 



-n 



(i,i,i) 



a,i,i) 



div,2 



(i,i,i) 



- * 



(i,i,i) 



(i,i,i) 



^(1,1,1) + ^(1,1,1) + ^(1,1,1) 



Now we can rewrite (16): 



U = EEE (%v,i,i,k*div,i,;,k + d div l2jl k*div, 2j ,k + ^n.,.k VI/ n.,k) 
where the divergence-free wavelets are simply obtained from the standard ones, for example: 



^(i,o,o) 
div,i 



j(i.o,o) 
div,2 



d (i,o,o) 



; (i,i,o)> 



rf (l,0,0) 



d 



(1,1,1) = 1. 

div,i 3 



(1,1,1) 



d 



(i,i,i) _ l 



-2d\ 



7(1.1,1) 



+ d 



(1,1,1) 



(i,i,i)^ 



rf (i,i,o) = d a,i,o) 

div,2 d 



div,2 



-d K 

3 \ «i 

The complement coefficients are in this case: 

_ j(1,0,0) \ , 1/^(1,0,0) 

MM U, 2MM-±M> ' 4^ 3MMM u 3MMM-i) 



( ,(1,0,0) _ ,(1,0,0) i^(i,o,o) _ ,(1,0,0) 

u n,fc ~~ a \MMM ' 4^ a 2MM '- n 



/(1,0,0) 



(dn) < 



d 



(1,1,0) _ Wj(1,1,0) 



(1,1,0^ 



n,fe — 2( i 

(i,i,i) _ i ^(i,i,i) +d (i,i,i) + d (i, 



2 

JU,1,1) _ 1/Jl' ^ 

^ u n ~~ 3 ^"2 "T "3 



+ ifd (1 ' 1,0) 
r 8V U 3,fci,fc 2 ,fe 3 



d (H,0) \ 



As for the two-dimensional case, the computation of divergence-free wavelet coefficients 
of any 3D vector field lies in a short linear combination of standard biorthogonal wavelet 
coefficients, arising from the fast wavelet transform. 
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3.2. Anisotropic divergence-free wavelet transforms 

In this section we will construct anisotropic wavelets that are divergence-free. Since the one- 
dimensional wavelets verify ipi = 4?/> , we derive easily divergence-free wavelet bases by 
tensor products of one-dimensional wavelets. We detail in the following the construction of 
such bases in the two and three dimensional cases. 



3.2.1. The anisotropic 2D case The 2D anisotropic divergence-free wavelets are given by: 

*df V j,k(^) 



2^ 1 {2^x l -k 1 )^(2^x 2 -k 2 ) 



where j = (ji, j 2 ) G Z is the scale parameter (which is different in both directions), and 
k = (fci, k 2 ) G Z 2 is the position parameter. When the indices k and j vary in Z 2 , the family 

^dfvj,k } forms a basis of Hdiv( M2 )- 
We introduced the complement functions: 



*Sj ik (*l,*2) 



2*^1(2*3:1 - fci)^o(2*x 2 - k 2 ) 
2^ {2^x 1 -k l )i, l {2^x 2 -k 2 ) 



The anisotropic divergence-free wavelet transform works similarly as the isotropic one 
but with fewer elements to be computed. The decomposition of a given vector function u 
begins with the anisotropic wavelet decomposition associated to the MRA (V^ 1 ® V®) x (V^° ® 
V--) (see section 2.3); 



-EE 

jeZ 2 keZ 2 



rf an ^an 



rf an ^an 



with: 



ij,k ij,k 2,j,k 2,j,k 



^i(2*xi- fciM)(2 J ' 2 x 2 - fc 2 ; 




* 2 j >k (^i^2) tf, (2* Xl - k l )^ 1 {2^x 2 - k 2 ) 

for j, k G Z 2 . Remark that for more simplicity, the dilated functions are not normalized, in 
L 2 -norm. 

u can be expanded onto the new basis: 

jeZ 2 keZ 2 
with the corresponding coefficients: 



d™ . 



2.12 



3 dTvj,k ~ 2^1+2^2 rf ij,k 2^1+2^2 a 2 j, k 

jan _ g£i jan , •m jan 
nj,k ~ 2^1+2^2 a ij,k + 2^1+2^2 a 2 ,j,k 



2-'i 



d™ 



(18) 
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3.2.2. The anisotropic 3D case In the same way, the (non normalized) anisotropic 3D 
divergence-free wavelets take the form: 



*dw,i,j,k^^) 



*dfv, 2 j,k(^,* 3 ) 



2^(2^ - h)M2 j2 X2 - Wo(2 J3 x 3 - fc 3 ) 
-2*Vo(2*zi - fci)Vi(2^ a x 2 - k 2 )MV 3 X3 ~ h) 




2^o(2 jl X! - fci)^i(2^ 2 x 2 - fc 2 )^ (2 J ' 3 x 3 - fc 3 ) 
-2^^o(2 il x 1 - h)M2 n X2 ~ Wi(2 3 ' 3 x 3 ~ h) 



V "dh,,j.k ( '-'^'- 



2'V;(2'\r, - fci)Vo(2 J ' 2 x 2 - A; 2 )^o(2 J3 X3 - fc 3 ) 



2^(2^ - Wo(2 J2 x 2 - fc 2 )^(2^x 3 - fc 3 ) 



with j = (ji, j 2 , j 3 ), k = (An, fc 2 , fca) e Z 3 . 

The 3D divergence-free basis is carried out by considering only two types of functions among 
the three above. As complement basis we introduce a function which is the most as possible 
orthogonal to the previous ones: 



2ii^ 1 (2*zi - A;i)^o(2 i2 x 2 - /c 2 )^ (2 J3 2 3 - k 3 ) 
2^ (2*xi - fc!)Vi(2^x 2 - fc 2 )Vo(2 J ' 3 x 3 ~ h) 
2^o(2 jl xi - h)M2 j2 X2 - k 3 )M^ 3 - k 3 ) 

The operations to compute divergence-free coefficients and complement coefficients are 
similar to the 2D case. 



4. An iterative algorithm to compute the Hodge wavelet decomposition 

4.1. Principle of the Hodge decomposition 

The Hodge decomposition consists in splitting a vector function u G (L 2 (R n )) n into its 
divergence-free component and a gradient vector. More precisely, there exist a pressure 
p and a stream-function ip such that: 

u = u^j v + Vp and u^j v = curl ip (19) 

Moreover, the functions curl tp and Vp are orthogonal in (L 2 (R n )) n . The stream-function ip 
and the pressure p are unique, up to an additive constant. 

In R 2 , the stream-function is a scalar valued function, whereas in R 3 it is a 3D vector function. 
This decomposition may be viewed as the following orthogonal space splitting: 

(L 2 (M"))" = H div (R")©H curl (M") 

where we note 

H div ( K ") = {v G (L 2 (M n ))7div v G L 2 (R n ), div v = 0} 
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the space of divergence-free vector functions, and 

H curl (K n ) = {v G (L 2 (M n ))7curl v G (L 2 (R n )) n , curl v = 0} 

the space of curl-free vector functions (if n = 2 we have to replace curl v G (L 2 (W n )) n by 
curl v G L 2 (M 2 ) in the definition). For the whole space R n , the proofs of the above decompo- 
sitions can be derived easily, by mean of the Fourier transform. In more general domains, we 
refer to [12]. 



The objective now is to provide with a wavelet Hodge decomposition. Since in 
the previous sections we have constructed wavelet bases of H ( jj v (E n ), we have to work 
analogously to carry out wavelet bases of H cur i (M n ). 



4.2. Construction of a gradient wavelet basis 

A definition of wavelet bases for the space H cur j (K n ) (n = 2, 3) has already been provided 
by K. Urban in the isotropic case [28]. We will focus here on the construction of anisotropic 
curl-free vector wavelets in the 2D case (it goes similarly in the n-dimensional case). 

This construction is very similar to the divergence-free wavelet construction, despite 
some crucial differences. The starting point here is to search wavelets in the MRA (Vj S*D 
Vj) x (Yj ® K/) mstea d of (Vj <g) Vj) x (Vj <g> V}), where the one-dimensional spaces V Q 
and V\ are related by differentiation and integration (proposition of section 2.4). 
Since H cur i(K 2 ) is the space of gradient functions in L 2 (lR 2 ), we construct gradient wavelets 
by taking the gradient of a 2D wavelet basis of the MRA iyj ® V^). If we avoid the L 2 - 
normalization, the anisotropic gradient wavelets are defined by: 



2^ (2^ Xl - £4)^1(2^x2 - k 
2^^1(2^x1 - k 1 )^ (2 j2 x 2 - k 



Thus,whenj = (j 1 ,j 2 )>k = (ki,k 2 ) vary in Z , the family {^^^j ^} forms a wavelet basis 
of H curl (R 2 ). 

The decomposition algorithm on curl-free wavelets {^^ r y j^} works similarly as the 
decomposition algorithm on anisotropic divergence-free wavelets. A vector function v is 
firstly approximated in a space (Vj <&Vj) x (Vj ® Vj ) by: 



1 = Ek e Z 2 c i,J,k ^o(2 J X! - k 1 )(j) l (2 J X2 - k 2 ) 
^keZ 2 c 2,J,k M 2Jx i ~ ki)(j)o{2 J X2 - k 2 ) 



By applying to vf the standard anisotropic wavelet transform of (Vj (g> Vj) and to vf the one 
of (V}®Vj), it rewrites: 



v*(x u x 2 ) 



v*{x u x 2 ) = £ Jlj2 <jE keZ 2 dJJ k M* 1 *! - fei) M^2 - k 2 ) 

v*(x u x 2 ) = E, 1J2 <,/E keZ 2 4j k M&xi - h) - h) 
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^ (2 J1 x 1 - fc 1 )^ 1 (2*z 2 - k 2 ) an ,# 

*2j,k^ 1 ' X2 J 





M2 n xi-ki)M2 j2 X2 



then 



J V iJ.k ij,k ^ 2j,k 2 j,k 

Thus, to compute the expansion of Pfv in terms of the gradient vector wavelets, we have to 
perform the change of basis: 




curlj.k 
^an 

ATJ,k 



2j!^,an,# _|_ 2j 2 \j/ an '# 
ij.k 2 j,k 
2i2v[/ an '# _ 2Ji\I/ an '^ 



ij,k 



2,j,k 



which leads to: 



(^curlj,k *curlj,k "Xj.k^ivj.k 



rf an :u* an : 



iU2<^keZ 2 

where the curl-free wavelet coefficients are obtained from the standard ones by: 

oh 

jan , z ran 



03i 9J2 

^an = z ^an , z jou 
curlj,k 2 2 J'i + 2 2 ^ Lj.k 2 2 ^ + 2 2 J2 2,j,k 



associated to complement coefficients: 

2 h 



d 



an 



d 



an 



2 ji 



7v,j,k 2 2 Ji + 2 2 ^ iJ.k 2 2 ^ + 2 2 ^ 2 J k 



(20) 



(21) 



(22) 



4.3. Implementation of the Hodge decomposition in the wavelet context 

From now on, our objective is to compute the wavelet decomposition of a given vector 
function v: this means to find a divergence-free component v^j v and an orthogonal curl-free 
component v cur | such that: 

v = v div + v curl 

where: 

v div = ^div.j.k^div.j.k v curl = ^ ^curl.j.k^curl.j.k 
j,k j,k 

are the wavelet expansions onto div-free and curl-free wavelet bases constructed previously 
(section 3.2.1 and 4.2). For more simplicity, we will focus on 2D anisotropic wavelet bases 
(and we will omit the superscript "an" in the notation of the basis functions). 

To provide with such decomposition, we have to overcome two problems: 
- The first one lies in the fact that div-free wavelets and curl-free wavelets form 
biorthogonal bases in their respective spaces, and as already noticed by K. Urban [28], they 
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would not give rise, in a simple way, to the orthogonal projections and v cm \ of v. As 
solution, we propose to construct, in wavelet spaces, two sequences ( v ^j v ) an d ( v c Ur i) mat 
will converge to and v cm \ ■ 

- The second difficulty is that div-free wavelets leave in spaces of the form (Vj <g> Vj) x 
(Vj ® Vj), whereas curl-free wavelets arise from (Vj (g) Vj) x (Vj (g) Vj), where V^ , Vq 
and Vq, Vq 1 are couples of spaces related by differentiation and integration. These spaces are 
different, and in order to construct our approximations ( v ^ v ) an d ( v ^ ur i)' we have to define 
a precise interpolation procedure between the two kinds of spaces. In particular, the spaces 
V^ , Vq can be suitably chosen from V^ , Vq. 

4.3.1. Iterative construction of the div-free and curl-free parts of a flow Let v = (vi,v 2 ) 
a vector function be given, and suppose that v is periodic in both directions, and known on 
2 J x 2 J grid points that are not necessarily the same for v\ and t> 2 . In the following, we will 
note: 

- Ijv an approximation of v in the space (Vj (g Vf) x (Vj ® Vj), given by some 
interpolating process. 

- 1* v an approximation of v in the space (Vj (g Vj) x (Vj (g Vj), also given by some 
interpolating process. 

We now define the sequences v^. y e (Vj <g Vj) x (Vj (g Vj) satisfying div v^. y = 0, and 

v curl G ^ X ® ^ satisf y in S curl v curl = °' as follows: 

- We begin with v° = Ijv and we compute v^- , the divergence-free wavelet 

decomposition of v°, and its complement v°, by formula (17,18): 

hv = V diy + V° = £ 4v j,k ^div,j,k + E < j,k *nj,k 

j,k j,k 
Then we compute at grid points the difference v — Vjj . 

Secondly we consider I*(v — v ^j v )> and we apply the curl-free wavelet decomposition 
(20,21,22), leading to a curl-free part and its complement: 

J *( v - V D = v curl + < = H rf curlj,k *curlj,k + £ ^.k^k 

j,k j,k 

Finally we define pointwise: v 1 = v — v^. y — v° uf j. 

- At step p, by knowing v p at grid points, we are able to construct a divergence free-part 
vJJ iy of Ijv p by (17), and v^, the curl-free component of if (v p - vj iy ) by (20) (v p - vj iy 
being computed at grid points). The next term of the sequence is again defined pointwise: 

v p+i = v p _ v p _ v p (23) 
div curl v ' 



Divergence-free Wavelets for Navier-Stokes 



22 



We iterate this process until ||v p ||£2 < e, and we obtain: 

p p 

curl 

p=i p=i 

= E(E4vj,k) %vj,k + E(E^urlj 1 k) *curlj,k 
j,k Vp =1 / j,k \p =1 f 

where the right hand side is an approximation of v, which interpolates the data up to an error 
e (e being given). 

For the moment, we are not able to prove theoretically the convergence to of 
the sequence (v p ): we will prove it experimentally in section 5.4, on arbitrary fields. 
Nevertheless, we can outline some remarks: 

- The convergence rate depends on the choice of complement functions ^ n j j^, ^ N j k- 

The more the L 2 -scalar products < ^iv j k' j' k' > anc * < ^curl j k' j' k' > are 
small, the faster the sequence converges. 

- Ideally, we would like to choose the interpolating operators Ij and I* such that the 
convergence doesn't depend on this choice. We propose below a choice for these operators, 
based on spline-quasi interpolation, which is satisfactory at relatively-slow convergence rate. 



4.3.2. Hodge-adapted spline interpolation In this part, we will detail our choice of operators 
Ij and I* , in the context of spline spaces of degree 1 (V®) and 2 (V^ 1 ) that we have introduced 
at the beginning. 

Let us suppose the components v\ and t> 2 of a velocity field v be known respectively at 
knot points 2~ J (A; 1 + |, /c 2 ) and 2~ J (ki, k 2 + h), for ki, fc 2 = 0.2 J — 1. This choice of grid is 
induced by the symmetry centers of scaling functions O of V and 4>i of Vi (see Figure 6). 




Figure 6. The two scaling functions of Vq and V\, and their symmetry centers 



For J given, Ij is chosen as an operator of quasi-interpolation (similarly to section 5.1.1) 
in the spline space (Vj <g> Vj) x (Vj® Vj) 
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^ v = E c k $ uk + E c k $ 2,j,k 
k k 

where $^ an d $2 are tne vec t° r scaling functions introduced in section 3.1.1. 

The second operator I* provides again with a quasi-interpolation of vector functions onto 
a new spline space (Vj ® Vj) x (Vj ® Vj ). Under interpolating considerations, we define: 

V° = {v ; v(x - 1/2) G Vb} = span{<p (x -1/2- k); k G Z} 

= {v ; v(x - 1/2) G Vi} = span{fa(x -1/2- k); k G Z} 
Hence we can write: 

k k 

where $1 jk and $2 jk 316 me anisotropic vector scaling functions of 3.1.1 built from 

0o = 0o and 0! = fa. 

5. Numerical experiments 

In this section, we present our numerical results concerning the application of divergence-free 
wavelet decomposition, for analyzing several data. We begin with the analyses of periodic, 
numerical, incompressible velocity fields in dimensions two and three, arising from pseudo- 
spectral codes. First, we have to take care of the initial interpolation of such fields, in 
order not to break the incompressible condition satisfied in Fourier space. Then, after the 
vizualisation of the divergence-free wavelet coefficients, we study the compression obtained 
through the wavelet decomposition. In the last part, we investigate and numerically prove 
the convergence of the algorithm presented in section 19, which provides with the wavelet 
Hodge decomposition of any flow. As an example, we compute the div-free component of the 
nonlinear term of the Navier-Stokes equations, and we extract the associated pressure, directly 
in wavelet space. In all the experiments, we will use divergence-free wavelets constructed with 
splines of degrees 1 and 2. 

5.1. Approximation of the velocity in spline spaces 

Usually, the data are provided by point values of the velocity field. The first step of the 
wavelet decomposition consists in interpolating the velocity coordinates on the suitable B- 
spline space. The arising problem is that this approximation may not conserve the divergence 
free condition that is verified in Fourier space, when velocity arise from a spectral code. 
The spline approximation of data, obtained through spectral methods, introduces a slight error 
for the divergence free condition. This difference may not be neglectable. For the turbulent 
fields we studied (2D and 3D) the error is about 1% of the L 2 -norm (i.e. 0.01 % of the energy). 

Thus we propose two ways to overcome the problem. The first way is to interpolate the 
velocity in the Fourier domain and to compute exactly its biorthogonal projection on wavelet 
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spaces. The second way is to interpolate on the divergence free B-spline spaces with a Hodge 
decomposition made by wavelet decompositions as it is proposed in the part 4.3.1 . This can 
be applied to any compressible flow. 

5.1.1. By quasi-interpolation The spline quasi-interpolation is a good compromise when 
we have to deal simultaneously with spline approximations of degree even and odd. In this 
context, the order of approximation is n + 1, by using B-splines of degree n [7]. An advantage 
of the procedure is that it may be applied in any case of boundary conditions. 

Let b be a B-spline scaling function (b = (p or <pi). Given the sampling f{k/N) 
(N = 2 J ), we want to compute scaling coefficients c*;, of a spline function f N , that will 
nearly interpolate the values f(k/N): 

f N (x) = Y,c k b(Nx-k) (24) 

fceZ 

/n is an interpolating function if: 

J2c k b(i-k) = f(^) WeZ 

fceZ 

For example, if we consider b = (pi (spline of degree 2), the previous condition implies: 

M^) = i(*-i + c«) = /(^) WgZ 

In order to avoid the inversion of a linear system, the quasi-interpolation introduces, instead 

of c e : 

c t = §[«£) + /(<±i)] - + /(<±?)] w 6 z 

By replacing q by q in (24), we obtain the following error at each grid point: 

i(c*_ 1 + ci) -/(I) = ![-/(? - 



2 V ^' J y N' 16 iV 



/ — 1 / / 4- 1 / 4- 2 

+ 4/(y)-6/(A) + 4/(i±l)-/(^)] 

f/ (4) W + 0(4) -with tfe] £ 



Therefore, the pointwise error of quasi-interpolation is order 4, for a sufficiently regular 
function. 



5.1.2. By using the Discrete Fourier Transform Since they are highly accurate, spectral 
methods are often considered as a reference technique for simulating incompressible turbulent 
flows. For periodic boundary conditions on the cube [0, l] 2 , the Discrete Fourier Transform is 
used to decompose the velocity u. 

If means the Discrete Fourier coefficients of u on a iV 2 regular grid, 

U k = Jp E u (]y) e N 

ne{o,i,...,N-i} 2 
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the velocity expansion in the Fourier exponential basis is: 

u(x)= ^k e2wk - X ( 25 ) 

ke{0,l,...,iV-l} 2 

In this context, the divergence-free condition div u = writes: 

k.u k = , Vk G {0, 1, . . . , iV — l} 2 (26) 

Assume now that the velocity field u we have to analyze (supposed to be 1 -periodic 
in both directions), is obtained from a spectral method and verifies the incompressibility 
condition in Fourier domain (26). To compute its decomposition in a divergence-free 
wavelet basis of M 2 , we have first to approximate u = (u\,U2) in the suitable space 
(Vj <g> Vj) x (Vj ® Vj) which has been introduced in section 3.1.1, where J corresponds to 
N = 2 J . Then we search for an approximate function uj = (uj i, uj 2 ) such that: 

v^2 J -l ^2 J -1 1 J. J. 

U J 1 ~~ 2^m=0 2^n 2 =0 C .J,n 1 ,n 2 rl, J.mVO, J,n 2 

2 J -1 v^2 J -l 2 



m=0 2^n 2 = 



C J,ni,n 2 V0,J,niVl,J,n 2 



For the choice of functions O an d 0i defined above (see equation (11)), the 
incompressibility condition div uj = takes the discrete form on the coefficients c l Jni m : 

C J,m,n 2 — C J,m+l,n 2 + C J,m,n 2 ~~ C J,m,n 2 +l = > ^( n l) n 2) (27) 

To conserve the incompressibility condition verified by u, a solution consists in 
considering uj as the biorthogonal projection onto the space (Vj (g> Vj) x (Vj 1 (g) Vj-), since 
we know that this projector commutes with the partial derivatives [19]. This is equivalent to 
consider that: 

C J,ni,n 2 = < u I 0*,J,m^O,J,n2 > 
C J,ni,n 2 = < U | 0o,J,ni^l,J,TO2 > 

Replacing u by its Fourier expansion (25), it follows: 

C \n u n 2 = ^k // 2 e 2mk - X ^ >J)ni (zi) 0S,J,n 2 (£2) cfeiGfe 2 



keio,!,---,^-!} 2 



kn 

Z27T- 



= 2" J ^ u k 0*(2- j 2ttA;i) 0*(2- J 2vrA;2) e 2 "'~ 
ke{o,i,...,iV-i} 2 

where §\ denote the (continuous) Fourier transforms of the dual scaling functions $>\, <$>\. 
Finally, we obtain an explicit form for the Discrete Fourier Transform DFT of the coefficients 
c J,n u n 2 (and in the same way for 4 ^): 



DFT(^ n ) k = u k 2-^(2^(27^)) 0S(2- J (27rfc 2 )) 

(28) 

DFT(c 2 >n ) k = u k 2- J ^(2- 7 (27r/c 1 )) ^(2- j (2tt/c 2 )) 
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It means that the discrete Fourier transform of coefficients c\ 



is given by the discrete 



Fourier transform of u, multiplied by tabulate values on [0, 2n] of the Fourier transform of the 
duals 02- In practice, we don't know the explicit forms of these functions, except by the 
infinite product: 



Nevertheless, the infinite product converges rapidly, which allows to obtain point values of 0q 
and </>*, with sufficiently accurate precision. 

In dimension 3, it goes similarly, by considering the biorthogonal projection of a 3D 
vector field u onto the space (Vj ® V] ® V]) x [V] <g> Vj ® V?) x (V$ <g> V] <g> V?). 

5.2. Analysis of 2D incompressible fields 

We focus in this part on the analyses of two-dimensional decaying turbulent flows. 

The first numerical experiment we present studies the merging of two same sign vortices. 
It concerns free decaying turbulence (no forcing term). The experiment was originally 
designed by M. Farge and N. Kevlahan [25], and often used to test new models [2, 13]. This 
experiment was here reproduced by using a pseudo-spectral/finite-difference method, solving 
the Navier-Stokes equations in velocity-pressure formulation. 

The initial state is displayed on figure 7 left. In a periodic box, three vortices with a 
gaussian vorticity profile are present; two are positive with the same intensity, one is negative 
with half the intensity of the others. The negative vortex is here to force the merging of the 
two positive ones. The time step was 5t = 1CT 2 and the viscosity v — 5 10~ 5 . The solution is 
computed ona512 x 512 point grid. 

The vorticity fields at times t = 0, t = 10, t — 20 and t = 40 are displayed on fig- 
ure 7. The last row of figure 7 displays the absolute values of the isotropic divergence-free 
wavelet coefficients at corresponding times, renormalized by 2 J at scale index j. As one can 
see, divergence-free wavelet coefficients concentrate on strong change in vorticity zones, that 
is around or in between vortices, or along vorticity filaments, as they are equivalent to second 
derivatives of the velocity. 

The second experiment deals with a decaying two-dimensional turbulent field, obtained 
with an initial state of random phase spectrum. That vorticity field was computed with the 
spectral code of [14], at a resolution 1024 x 1024, and a Reynolds number of 3.5 x 10 4 . This 
field has been kindly provided to us by G. Lapeyre [18]. Figure 8 left represents the vorticity, 



k (0 (^) = 



0o(On,>i(2-cos(^)) 




n j > 1 (2 - cos(£2- J )) 
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after 40 turnover time-scale, where it exhibits the emergence of coherent structures together 
with strong filamentation of the flow field outside the vortices. 

We show in Figure 8 right, the isotropic divergence-free wavelet coefficients of the 
corresponding velocity field, in L inf -norm. 

As expected, the wavelet coefficients get an insight into the energy distribution over 
the scales of the flow. As one can see on Figure 8, the energy at smallest scale (or highest 
wavenumbers) is localized along the strong deformation lines, and fits the filamentation 
between vortices, or with strong changes in vortices. The top-right square corresponding 
to vertical isotropic wavelets O&^J, ^ ) exhibits vertical structures, whether the bottom-left 

square corresponding to horizontal wavelets (^.'^ . , ) exhibits horizontal deformation lines. 



Now we investigate the compression properties of the divergence-free wavelet analysis: 
as predicted by the nonlinear approximation theory (see [3]), the compression ratio in energy- 
norm is governed by the underlying regularity of the velocity field in some Besov space. 
Let u an incompressible field be given, its divergence-free wavelet expansion writes: 

u — u + (^div,j,k ^div,j,k ^div,i,k ^div,j,k ^divj.k ^div,j,k) 

^° k G Z 2 

The nonlinear approximation of u relies on computing the best N-terms wavelet 
approximation by reordering the wavelet coefficients: 

\d% i I > l<r i I > • • • > W A N - 1 I > ... 
1 divji.ki 1 1 divj 2 ,k 2 ' 1 divjjv.kjv 1 




Figure 7. Vorticity fields at times t = 0, t = 10, t = 20 and t = 40, and corresponding 
divergence-free wavelet coefficients of the velocity. 
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and introducing 

N 

^(u) = Uo + ^^ V) . ki ^ v . ijki (29) 

i=i 

Then we have 

\\u-E N (u)\\ L 2<c(^] Hull^ (30) 



iV 



</ 



if the quantity ||u||^ s , 9 = X] e jk ^div k * s ^ n ^ te ' w ^ tn \ = \ + n means t na t u 
belongs to the Besov space B s q ,q ). As stated in [3], the evaluated regularity s can't be larger 
than the order of polynomial reproduction in scaling spaces plus one (that equals the number 
of zero moments of the dual wavelet). In our experiment, the dual spline wavelets i/j^ and 
ijj* introduced in 2.2 have respectively two and three zero moments, which only allows us to 
evaluate regularities smaller than two. 

Figure 9 shows the nonlinear compression of divergence-free wavelets, provided on the 1024 2 
turbulent field. The curve represents the L 2 -error ||u — £jv(u)|| L 2, versus N, in log-log plot. 
The convergence rate measured on the curve is s ~ 1.35 , which induces that the velocity flow 
belongs to the corresponding Besov space B s q ,q with q = 0.85. 
When looking at the compression curve on figure 9, we observe three zones: 

- First, large scale wavelets capture the large scale structure of the flows. Consequently, 
the compression progresses slowly and irregularly. 

- Then we observe a linear slope that represents the nonlinear structure of the turbulent 
flows. In this region, we are able to evaluate the regularity of the field. 

- The last region corresponds to an abrupt decrease, due to the fact that the data are 
discrete. 




Figure 8. Vorticity field for a 1024 x 1024 simulation of decaying turbulence (left), and the 
corresponding divergence-free wavelet coefficients of the velocity field (right). 
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One can also remark on Figure 9 that only 1.2% of the coefficients recover about 99% of 
the L 2 -norm. 




Figure 9. L 2 -error provided by the nonlinear N-best terms wavelet approximation (29): in 
log-log plot, L 2 -error (30) versus N for a 2D turbulent flow. 

The same experiment was carried out on the three interacting vortices, but due to the few 
number of vanishing moments of the wavelets we use (2), the slope of the curve saturates at 
s — 2, meaning that these fields are more regular. 

5.3. Analysis of a 3D incompressible field 

In this part we consider a three dimensional periodic field, arising from a freely decaying 
isotropic turbulence, and kindly provided to us by G.-H. Cottet and B. Michaux [4]. 
The experiment deals with an initial velocity condition of Gaussian distribution, and 128 3 
collocation points. Figure 10 displays the vorticity isosurfaces corresponding to about 40% of 
the maximum vorticity at five turnover times. 




1G '1C 



Figure 10. Isosurface of vorticity magnitude after 5 large-eddy turnovers provided by a 
spectral method [4] . 

The divergence-free wavelet decomposition of the corresponding velocity field is com- 
puted, and displayed on Figures 11 and 12. As explained in section 3.1.2, the isotropic 
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3D divergence-free wavelet decomposition provides with 14 generating wavelets 1 
div,2j,k 
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Figure 1 1 left shows the corresponding renormalized coefficients 2 J c/ ( jj v 1 j^, whereas 

Figure 11 right shows the 2- 7 ^ V2j .] { , until j = 6. The smallest scale (j = 7) wavelet 

coefficients are displayed on Figure 12 below, for two kinds of generating wavelets: we choose 

which corresponds to horizontal structures, and vj/^ .' ' 1 ^ which exhibits vertical ones. 
div,2 r div,i 




Figure 11. Isosurface 0.2 of divergence-free wavelet coefficients associated to 1 . ^ (left) 
and to 2 . k (right), in absolute value. 

Figure 13 displays the nonlinear compression error: we have computed the convergence 
rate on the linear part of the graph (which is shorter by comparison with the 2D case, due the 
low resolution) and we have found s « 1.45. 

5.4. Analysis of 2D compressible fields 

We presented in section 4.3, an algorithm which gave rise to a wavelet Hodge decomposition 
of any flow. In order to numerically prove that it always converges, we have tested the method 
on various random two-dimensional fields. We constructed some of them by summing random 
gaussians, and we modified their Fourier spectra in order to vary the regularity. Figure 14 
displays the L 2 -norm of the residual, in terms of the number of iterations, for four different 
vector functions. We can infer the following conclusions: 

- For all functions we have tested, the method converges, and curves shows that, except 
at the early beginning, the convergence is exponential. 

- The slope of the curve do not depend to much on the number of grid-points, but the 
curve itself, corresponding to 1024 2 grid points, is upon these of 256 2 grid points. 
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(1,0,0) 

div,2,j,k 



Figure 13. L 2 -error provided by the nonlinear N-best terms wavelet approximation (29): in 
log-log plot, L 2 -error (30) versus N for a 3D turbulent flow. 

- The convergence rate increase with the number of vanishing moments of the dual 
wavelets. 

In the futur, we will investigate the influence of the wavelet bases and of the interpolating 
projectors, on the convergence rate. 

Since our main objective for further research is to use divergence-free wavelets 
for solving the Navier-Stokes equations, we have to provide with the wavelet Hodge 
decomposition of the nonlinear term (u.V)u. Indeed, although u should be incompressible, 
this term breaks the divergence-free condition and yields a compressible part. As an 
illustration, we consider as u the 2D turbulent field displayed on Figure 8, and we compute, 
by mean of our wavelet Hodge decomposition, the div-free and curl-free wavelet components 
of (u.V)u. Figure 15 shows the anisotropic wavelet coefficients of the divergence-free part 
(left) and of the curl-free part (right) of the (u.V)u arising from this decomposition. 

Figure 16 (left) displays the vorticity field associated to the divergence-free part of 
(u.V)u, while Figure 16 (right) represents the pressure issued from the curl-free term, that is 
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256x256 points 
2 zero moments 



256x256 points 
more regular 
2 zero moments 



1024x1024 point 
2 zero moments 



256x256 points 
3 zero moments 



Number of iterations 




Figure 15. Anisotropic wavelet coefficients corresponding to the wavelet Hodge 
decomposition of (u.V)u: divergence-free coefficients (left), and curl-free coefficients (right). 



easily reconstructed in wavelet domain, as it will be explain below. 

The extracted div-free part of (u.V)u is all what is needed to compute the time-evolution 
of the velocity in the incompressible Navier-Stokes equations. Thanks to the curl-free wavelet 
definition, we are also able to directly reconstruct the pressure from the curl-free coefficients 
of Vp: 



Divergence-free Wavelets for Navier-Stokes 



34 




Figure 16. Vorticity (on the left) and pressure (on the right) derived from the wavelet Hodge 
decomposition of the nonlinear term uVu, with u displayed on Figure 8. 



Indeed, with periodic boundary conditions, the curl-free part of (u.V)u writes: 

J-l 2-n-l 2*2-1 



[(u.V)u] curl 1 



dp 
dx 



ji,32=0 fei=0 fc2=0 
J-l 2*1-1 



ji=0 fei=0 



[(u.V)u] curl 



dp 
dy 



J-l 2*1-12*2-1 

E E E <Wi,* 2 2J2 ^ ^ - fci) ^b(2* V - h) 

J1J2=0 fcl=0 fc2=0 

J-l 2*2-1 

+ E E <SuiU 2" *>(2* y - A*) 

j 2 =0 fc 2 =0 

Be integrating the system (we recall that r tp[ = 4 ip , see the definition of gradient 
wavelets), we obtain (up to a constant): 

J-l 2*1-12*2-1 

4 p(s, y) = E E E d cur\M,kM 2jl x - fc ^C 2 * f - **) 
iu'2=o fci=o fc 2 =o 

+ EE rf curU>(2 j1 --^ 

ji=0 fci=0 
J-l 2*2-1 

+ EE d luA M ^ h v-^) 

j 2 =0 fc 2 =0 
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Thus the computation of the pressure is no more than a standard anisotropic wavelet 
reconstruction in Vj x Vj, from the curl-free coefficients. By comparison to the pressure 
computed in Fourier domain, we have found a relative error of 2.5 10~ 4 in the L 2 -norm, which 
probably arises from the interpolating process. On the other hand the difference between the 
Leray projection (in Fourier space) and the wavelet projection onto the divergence-free space 
represents 1% of the L 2 -norm, that is to say 0.01% of the energy. Figure 17 displays the 
localisation of this error. 



y 



v 



v. 



:\ v 



Figure 17. Error between the divergence-free part of (u.V)u (obtained through Fourier 
transform) and the one provided by the wavelet Hodge decomposition. 



Conclusion and perspectives 

We have presented in detail the construction of 2D and 3D divergence-free wavelet bases, and 
a practical way to compute the associated coefficients. We have introduced anisotropic div- 
free and curl-free wavelet bases, which are more easy to handle. We have shown that these 
bases make possible an iterative algorithm to compute the wavelet Hodge decomposition of 
any flow. Thus, numerical tests prove the feasibility of divergence-free wavelets for simulating 
turbulent flows in two and three dimensions. A divergence-free wavelet based solver for 2D 
Navier-Stokes equations is underway and will be reported in a forthcoming paper. 

An important issue that must be addressed is the great ability of the method: although 
all numerical tests have been presented in the periodic case, the method extends readily to 
non-periodic problems, by using wavelets adapted to the proper boundary conditions [28, 24], 
in the div-free construction. Another point is since we consider the (u, p) -formulation for the 
Navier-Stokes equation, and since we are able to compute the Leray projector in the wavelet 
domain, the method extends easily to the 3D case. At last, this method should be competitive 
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by comparison to a classical Fourier method in the non-periodic case: indeed, the periodic 
case corresponds to boundary conditions for which spectral methods are obviously fast, while 
it is clear that wavelet methods take advantage both of the compression properties of the 
wavelet bases for functions and for operators, in any case. 

Acknowledgments 

The authors would like to thank G.H. Cottet and G. Lapeyre for helpfully providing to them 
numerical turbulent flows for analyses. This work has been supported in part by the European 
Community's Human Potential Programme under contract HPRN-CT-2002-00286, "Breaking 
Complexity". 

References 

[1] C.-M. Albukrek, K. Urban, W. Dahmen, D. Rempfer, and J.-L. Lumley, Divergence-Free Wavelet 

Analysis of Turbulent Flows, J. of Scientific Computing 17(1): 49-66, 2002. 
[2] P. Charton, V. Perrier, A pseudo-wavelet scheme for the two-dimensional Navier-Stokes equations, 

Comp. Appl. Math. 15(2): 139-160, 1996. 
[3] A.Cohen, Wavelet methods in numerical analysis, Handbook of Numerical Analysis, vol. VII, 

P.G.Ciarlet and J.L.Lions eds., Elsevier, Amsterdam, 2000. 
[4] G.-H. Cottet, B. Michaux, S.Ossia and G. Vanderlinden, A comparison of spectral and vortex 

methods in three-dimensional incompressible flows, J. Comp. Phys, 175, 2002. 
[5] W. Dahmen, A. Kunoth and K. Urban, A wavelet-Galerkin method for the Stokes problem, 

Computing 56, 1996, 259-302. 
[6] I. Daubechies, Ten lectures on Wavelets, SIAM book, Philadelphia, Pennsylvania, 1992. 
[7] C. De Boor, A Practical Guide to Splines, book, Springer- Verlag New York Inc., 2001. 
[8] M. Farge, Wavelet transforms and their applications to turbulence, Ann. Rev. Flu. Mech. :395- 

457, 1992. 

[9] M. Farge, N. Kevlahan, V. Perrier & E. Goirand, Wavelets and turbulence, Proc. IEEE 84(4), 
639-669,1996. 

[10] M. Farge and K. Schneider, Coherent Vortex Simulation (CVS), A Semi-Deterministic Turbulence 

Model Using Wavelets, Flow, Turbulence and Combution, 66: 393-426, 2001. 
[11] J. Frohlich and K. Schneider, Numerical simulation of decaying turbulence in an adaptive 

wavelet basis, Appl. Comput. Harmon. Anal., 3: 393-397, 1996. 
[12] V. Girault, PA. Raviart, Finite element approximations of the Navier-Stokes equations, Lecture 

Notes in Mathematics, Springer- Verlag, 1979. 
[13] M. Griebel and F. Koster, Adaptive wavelet solvers for the unsteady incompressible Navier-Stokes 

equations, Advances in Mathematical Fluid Mechanics, J. Malek and J. Necas and M. Rokyta 

eds, Springer- Verlag, 2000. 
[14] B.L. Hua and D. Haidvogel, Numerical simulations of the vertical structure of quasi- geostrophic 

turbulence, J. Atmos. Sci. 3:2923-2936, 1986. 
[15] J. Ko, A.J. Kurdila and O.K. Rediniotis, divergence-free Bases and Multiresolution Methods for 

Reduced-Order Flow Modeling, AIAA Journal, 38(2): 2219-2232, 2000. 
[16] F. Koster and M. Griebel and N. Kevlahan and M. Farge and K. Schneider, Towards an adaptive 

wavelet-based 3D Navier-Stokes solver, Numerical flow simulation I, Notes on Numerical Fluid 

Mechanics, Vol. 66, 339-364, E.H. Hirschel eds, Vieweg-Verlag, Braunschweig, 1998. 



Divergence-free Wavelets for Navier-Stokes 



37 



[17] J.-P. Kahane and P.-G. Lemarie-Rieusset, Fourier series and wavelets, book, Gordon & Breach, 
London, 1995. 

[18] G. Lapeyre, Topologie de melange dans un fluide turbulent geophysique (in french), These de 

doctorat de l'Universite Paris VI, 2000. 
[19] P.-G. Lemarie-Rieusset, Analyses multi-resolutions non orthogonales, commutation entre 

projecteurs et derivation et ondelettes vecteurs a divergence nulle (in french), Revista 

Matematica Iberoamericana, 8(2): 221-236, 1992. 
[20] P.- G. Lemarie-Rieusset, Un theoreme d'inexistance pour les ondelettes vecteurs a divergence 

nulle (in french), C. R. Acad. Sci. Paris, t. 319, Serie I, p. 811-813, 1994. 
[21] J. Lewalle, Wavelet transform of the Navier-Stokes equations and the generalized dimensions of 

turbulence, Appl. Sci. Res. 51(1-2): 109-1 13, 1993. 
[22] S. Mallat, A Wavelet Tour of Signal Processing, book, Academic Press, 1999. 
[23] C. Meneveau, Analysis of turbulence in the orthonormal wavelet representation, Journal of Fluid 

Mechanics 232: 469-520, 1991. 
[24] P. Monasse and V. Perrier, Orthonormal wavelet bases adapted for partial differential equations 

with boundary conditions, SIAM J. on Math. Analysis 29(4): 1040- 1065, 1998. 
[25] K. Schneider, N. Kevlahan, M. Farge, Comparison of an adaptive wavelet method and nonlinearly 

filtered pseudo-spectral methods for two-dimensional turbulence, Theor. Comput. Fluid Dyn. 

9: 191-206, 1997. 

[26] K. Urban, A Wavelet-Galerkin Algorithm for the Driven— Cavity— Stokes— Problem in Two Space 

Dimensions, RWTH Aachen, Preprint 1994. 
[27] K. Urban, Using divergence free wavelets for the numerical solution of the Stokes problem, 

AMLF96: Proceedings of the Conference on Algebraic Multilevel Iteration Methods with 

Applications, 2: 261-277, University of Nijmegan, The Netherlands, 1996. 
[28] K. Urban, Wavelet Bases in H(div) and H(curl), Mathematics of Computation 70(234): 739-766, 

2000. 

[29] K. Urban, Wavelets in Numerical Simulation, Springer, 2002. 



Divergence-free Wavelets for Navier- 
Figure captions 



Divergence-free Wavelets for Navier-Stokes 



39 



Figure 1. From left to right: the scaling function 4> with its associated symmetric wavelet with 
shortest support, and their duals: the dual scaling function <fi* and the dual wavelet ip*. 



Figure 2. Scaling functions and associated wavelets with shortest support, for splines of 
degree 1 (left) and 2 (right). 



Figure 3. Anisotropic 2D wavelet transform. 



Figure 4. Isotropic 2D wavelet transform. 



Figure 5. Isotropic 2D generating divergence free wavelets (left), ^jj 1 ' (center) and 

(right). 



Figure 6. The two scaling functions <pQ and <pi, and their symmetry centers. 



Figure 7. Vorticity fields at times t = 0, t = 10, t = 20 and t = 40, and corresponding 
divergence-free wavelet coefficients of the velocity. 



Figure 8. Vorticity field for a 1024 x 1024 simulation of decaying turbulence (left), and the 
corresponding divergence-free wavelet coefficients of the velocity field (right). 
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Figure 9. L 2 -error provided by the nonlinear N-best terms wavelet approximation (29): in 
log-log plot, L 2 -error (30) versus N for a 2D turbulent flow. 



Figure 10. Isosurface of vorticity magnitude after 5 large-eddy turnovers provided by a 
spectral method [4]. 



Figure 11. Isosurface 0.2 of divergence-free wavelet coefficients associated to 1 . ^ (left) 
and to ^^jy 2 . k (right) in absolute value. 



Figure 12. Isosurface 0.06 of divergence-free wavelet coefficients associated to vf', 1 / ' '' 
& & div,2,j,k 

(left) and to ^^jy'i^jk: ^^ht)> m absolute value. 



Figure 13. L 2 -error provided by the nonlinear N-best terms of wavelet approximation (29): in 
log-log plot, L 2 -error (30) versus N for a 3D turbulent flow. 



Figure 14. Convergence curves of the iterative wavelet Hodge algorithm. 



Figure 15. Anisotropic wavelet coefficients corresponding to the wavelet Hodge 
decomposition of (u.V)u: divergence-free coefficients (left), and curl-free coefficients (right). 



Figure 16. Vorticity (on the left) and pressure (on the right) derived from the wavelet Hodge 
decomposition of the nonlinear term uVu, with u displayed on Figure 8. 



Figure 17. Error between the divergence-free part of (u.V)u (obtained through Fourier 
transform) and the one provided by the wavelet Hodge decomposition. 




number of coefficients 



